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Abstract. Recent experimental observations of time-dependent beatings in the two- 
dimensional echo-spectra of light-harvesting complexes at ambient temperatures have 
opened up the question whether coherence and wave-like behaviour plays a significant 
role in photosynthesis. We perform a numerical study of the absorption and echo- 
spectra of the Fenna-Matthews-Olson (FMO) complex in chlorobium tepidum and 
analyse the requirements in the theoretical model needed to reproduce beatings in 
the calculated spectra. The energy transfer in the FMO pigment-protein complex is 
theoretically described by an exciton Hamiltonian coupled to a phonon bath which 
account for the pigments electronic and vibrational excitations respectively. We use 
the hierarchical equations of motions method to treat the strong couplings in a non- 
perturbative way. We show that the oscillations in the two-dimensional echo-spectra 
persist in the presence of thermal noise and static disorder. 
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1. Introduction 

Light liarvesting complexes are the part of photosynthetic systems that channel energy 
from the antenna to the reaction centre. One of the most studied complexes is the 
Fenna-Matthews-Olson (FMO) complex [H El |3] which is part of the photosynthetic 
apparatus of green sulphur bacteria. The observation of long lasting beatings, up to 
1.8 ps at T = 77 K in time-resolved two-dimensional (2d) spectra [U El E] has generated 
enormous interest in theoretical modelling the process of energy transfer through light- 
harvesting complexes. Recently molecular-dynamics simulations started to develop a 
more microscopic picture of the dynamics [3 |8], but an atomistic and fully quantum- 
mechanical model of the whole process reproducing the experimentally obtained spectra 
is not available. Therefore theoretical models treat the light-harvesting pigment-protein 
complex as an exciton system coupled to a bath which accounts for the vibrational 
modes of the pigments. 

The choreography of the exciton-energy transfer can be elucidated with 2d echo- 
spectroscopy using three laser pulses which hit the sample within several femtoseconds 
time spacings. The experimentally measured 2d echo-spectra of the complex show 
a variety of beating patterns with oscillation periods roughly consistent with the 
differences in eigenenergies of the excited system obtained from fits to the absorption 
and 2d echo-spectra [5]. At long delay times between pump and probe pulses the system 
displays relaxation to the energetically lowest lying exciton eigenstates as expected from 
the interactions between the electronic degrees of freedom and the vibrational modes of 
the complex which drive the system into a state of thermal equilibrium [9] . 

One main goal of the research in this field is to clarify the role of a disordered 
and fluctuating environment in the energy-transfer efficiency. Theoretical quantum- 
transport studies show that the energy transfer in these systems results from a balance 
between coherent and dissipative dynamics [TOl|lT]. Calculations at room temperature 
also show that the optimal energy-transport efficiencies appear as a compromise between 
coherent dynamics and thermal fluctuations [12]. 

Another open question in the theoretical study of these systems is the description 
of experimental observables, i.e. the absorption and the 2d echo-spectra. Most of the 
theoretical studies analyse the beatings in the population dynamics and compare them 
with the ones observed in the experimentally measured 2d echo-spectra. However, the 
2d echo-spectrum reflects different information than the population dynamics and needs 
to be calculated and analysed separately in detail. In principle techniques such as 
quantum state and process tomography made out of a sequence of 2d echo-spectra 
can be used to map out the complete density matrix ^3]. In contrast to energy- 
transfer efficiency studies where an initial excitation enters the complex at specific sites 
close to the antenna, in 2d echo-spectra the whole complex is simultaneously excited. 
Also the two-exciton manifold yields prominent contributions to the signal, resulting 
in negative regions in the 2d echo-spectra. The non-pertubative calculation of 2d 
echo-spectra presents a considerable computational challenge due to the presence of 
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two excitons giving raise to excited state absorption and the requirement to consider 
ensemble averages over differently orientated complexes with slightly varying energy 
levels. Previous calculations have employed Markovian approximations [T5] or 
exclude the double exciton manifold [TB] . In addition the systematic study of beatings 
in a series of 2d echo-spectra requires effective means to calculate a huge number 
of such spectra. So far no theoretical method has been able to describe the long- 
lasting beatings in the time- resolved 2d spectra [T71 HH |TB] . One possible explanation 
for the persistence of long coherence times has been the sluggish absorption of the 
reorganisation energy by the molecule, which requires theoretical descriptions that go 
beyond the Markovian approximation and the rotating wave approximation [18]. The 
hierarchical equations of motions (HEOM), first developed by Tanimura and Kubo [IS] 
and subsequently refined in [201 ED [22l |23], show oscillations in the dynamics of the 
exciton populations that persist even at temperature T = 300 K [2^ 125] . The HEOM 
include the reorganisation process in a transparent way and are directly applicable for 
computations at physiological temperatures. A recent calculation at temperature 77 K 
of the 2d echo-spectra with the HEOM method has been performed by Chen et al [T7] , 
which does not display the strong beatings observed experimentally. Chen et al conclude 
that the agreement between theory and experiments needs to be improved. Here, 
we analyse under what conditions the theoretical calculations give better agreement 
with the experimental results. We calculate the absorption and 2d echo-spectra in a 
temperature range between 77 K to 277 K using the HEOM and investigate the role of 
pigment and protein vibrations in the calculated spectra. We systematically discuss 
the parameters needed to obtain results that agree better with the experimentally 
measured spectra. We discuss the role of static disorder for the shapes of the peaks and 
the appearance of beatings in the time-resolved 2d echo-spectrum and their robustness 
against temperature changes and static disorder. 

The paper is organised as follows: the exciton model and the HEOM are introduced 
in sec. [2J The absorption spectra and the line shapes obtained from the HEOM model 
of the FMO complex are analysed in sec. [3] where we compare different methods for 
performing rotational averages. In sec. HI we discuss the 2d spectrum and analyse how 
disorder affects it. Furthermore we analyse the persistence of beatings in the spectra 
as a function of the delay time for several temperatures and with and without static 
disorder. 

2. Exciton model of the FMO complex 

The FMO complex described by Fenna, Matthews and Olson [21 |3] is a pigment- 
protein complex that channels the energy in the photosynthetic apparatus in green 
sulphur bacteria. The FMO complex is arranged as a trimer with the different subunits 
interacting weakly with each other. We thus restrict our study to a single subunit which 
contains eight bacteriochlorophyll molecules (BChls) [2^1 EZ]- The BChls pigments, 
which are wrapped in a protein environment, can be electronically excited and are 
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coupled due to dipole-dipole interactions. The BChls form a network that guides 
the energy from the antenna to the reaction centre. Since the eighth BChl is only 
loosely bound it usually detaches from the others when the system is isolated from its 
environment to perform experiments [SHIES]. Therefore, we do not take the eighth BChl 
into account and describe the energy transfer using the seven site Frenkel exciton model 
[29| l30l [3T| . The single exciton manifold is given by 

Hl-=f^e,\k){k\ + J2j,d\k){l\ + \l){k\) (1) 

k=l k>l 

with A'^sites = 7. The state \k) corresponds to an electronic excitation of BChlfc 
with site energy e^. Dipole-dipole couplings between different chromophores yield the 
inter-site couplings Jki, which lead to a delocalisation of the exciton over all BChls. 
The site energies of the BChls = e° + consist of 'zero-phonon energies' e°, 
and a term due to the reorganisation Hamiltonian ifreorg = '^k=T ^k\k) {k\. The 
reorganisation energy Afc = /'^)k denotes the energy which is removed 

from the system to adjust the vibrational coordinates to the new equilibrium value 
[3H [211 [32] . Here, denotes a dimensionless displacement of the ith oscillator measured 
in units of \Jh/{rniUi), where mj and Ui are mass and frequency of the ith oscillator. 
The vibrational modes of the BChl are modelled as a set of harmonic oscillators 
-f^phon = ^ihujih\hi. The vibrational environment couples linearly to the exciton 
system ifex-phon = Yl!k=T + b])^ \k) {k\ and induces fluctuations in the 

site energies. We neglect correlations between the vibrations at different sites and model 
the coupling strength by a structureless spectral density 

A^) = 2A^^, (2) 

where we assume identical couplings = A at each site. This choice of spectral 
density facilitates the use of the HEOM method since the time-dependent bath 
correlations then assume an exponential form. For the FMO complex a variety of 
spectral densities are discussed in the literature, where some are a parametrisation 
of experimental fluorescence spectra [33], while others are extracted from molecular 
dynamics simulations [M] . We use the model Hamiltonian for chlorobium tepidum given 
in Ref. [T^. The Hamiltonian originates from previous calculations in |35l[36]- Cho et 
al. [2] changed the coupling constant between BChl 5 and BChl 6 and obtained new site 
energies from a fit to experimental results. The Hamiltonian of Ref. [H] has been used 
in previous studies of 2d echo-spectra [HI [161 IH] ^-iid facilitates comparisons between 
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the different approaches. In the site basis {\k)} the matrix elements read 
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in units of cm^ , where we have added an energy shift As = 12075 cm^ to the diagonal 
elements of the exciton Hamiltonian H^'^^ to shift the positions of the eigenenergies to 
the frequency range of the figures in Ref. [H]. We set the parameters of the spectral 
density to A = 35 cm^^ and 7^^ = 50 fs. These parameters are chosen to provide 
best agreement to the smooth part of the spectral density in Ref. [33], where the low 
frequency component of the spectral density was fitted to fluorescence line narrowing 
spectra ^7\. Since a single Lorentz-Drude peak is not sufficient to accurately reproduce 
the smooth part of the spectral density of Ref. [33] , we choose the parameters to obtain 
agreement within a factor of 1.2 to 1.5 in the energy range from 100 to 150 cm~^, which 
corresponds to the typical difference between the exciton energies. 

The HEOM approach rewrites time non-local effects into a hierarchy of coupled 
equations of motion for a set of auxiliary matrices a". The superscript n is a vector 
of Ai'sites entries, where gives the order to which the bath of BChU is taken into 
account. The reduced density matrix p{t) = trphonon -R(^) describing the exciton system 
is obtained after tracing out the vibrational degrees of freedom. The time evolution of 
the reduced density matrix is given by 



dt 
with 

d . 



p (t) = -i£eP (t) + J2 ('^) 



(4) 
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i£e +5^7^^) (t) 

k=l 



J2 ^kCT'''"''' (t) + ^kQka''-'" it) . (5) 

k=l k=l 

Ce is the Liouville operator of the coherent exciton dynamics and we define Cfc as the 
kth unit vector with iVsitcs = 7 entries. The operators and 6fc are defined by their 
action on a test operator A 



^kA 
QkA 



i^dfc) {k\A + A\k) {k\) 



(6) 
(7) 



where A 



[\k){k\^ A\. The auxiliary matrices are initialised to a"^ (0) = for 7^ 
and (T° (t) = p (t) ^24j. The hierarchy is truncated for sufficiently large A^max = J2k=T '^k, 
where the diagonal coupling in ([5]) becomes dominant. For the parameters used in 
the calculation, reorganisation energy A = 35 cm~^ and phonon relaxation time scale 
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7^^ = 50 fs, we verified that the relative differences in the interesting energy region 
of the absorption spectrum using A^max = 4 versus A^max = 12 are below 0.13 percent. 
In the 2d echo-spectrum we compared the rephasing stimulated emission diagram for 
A'max = 4 and iV^ax = 12 and found a difference of less than 0.07 percent of the 
amplitude. We therefore consider a truncation at NrascK = 4 sufficient to reach converged 
results. The derivation of the hierarchy stated in (jl]) and ([5]) relies on a high temperature 
approximation. For a temperature of 77 K the high temperature approximation is not 
valid and low temperature correction terms have to be included [25]. Within the low 
temperature correction we replace 

^ sites o \ O 

2^fe 27 ^^x^^x 

where vi = 2tt/ (ih is the first bosonic Matsubara frequency following [25]. Agreement 
with [5H] and proper thermalization to the Boltzmann distribution [12] confirms that 
the first order correction in ([8]) is indeed sufficient. A further improvement of the low 
temperature correction can be attained by including higher Matsubara terms. Recently, 
also a Fade spectrum decomposition has been proposed [391 SO] which shows faster 
convergence than the Matsubara decomposition. 

The numerical evaluation of the HEOM requires to propagate a large number of 
auxiliary matrices. Although the number of auxiliary matrices can be reduced by a 
scaled HEOM approach [1T| WI\ as well as by advanced filtering techniques [H], the 
HEOM are computationally very demanding, which sets limits to treatable system 
sizes. For 2d echo-spectroscopy one has to propagate the set of auxiliary matrices 
not only for a huge number of time steps (typically of the order of 10^), but also the 
dimension of the auxiliary matrices gets enlarged since excited state absorption requires 
to extend the exciton Hamiltonian to the two exciton manifold. We overcome the 
computational limitation by using a GFU (Graphics Frocessing Unit) implementation 
where the computation time goes down with the number of cores on a single GFU ^12j . 
For example on a single NVIDIA C2050 graphics board we obtain a 450 fold speedup 
compared to an equivalent single-core GFU implementation. With our efficient GFU 
algorithm the computation time for a 2d spectrum is reduced to 50 minutes on a single 
GFU. 



3. Absorption spectra of the FMO complex 

Spectroscopy provides a tool that gives direct insight into the energy states of a quantum 
system. Absorption spectroscopy is used to identify the energy eigenstates and their 
interaction strengths with the radiation field, expressed in terms of dipole moments. 
The energy eigenstates and the dipole moments are commonly identified with the peak 
position and the peak areas respectively. The spectral peaks are broadened due to 
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the interaction of the energy eigenstates with the remaining degrees of freedom of the 
system. For an ensemble of hght-harvesting complexes the coupling of the electronic 
dynamics to the vibrational modes and the presence of fluctuations in the energy levels, 
due to e.g. the motion of the surrounding proteins, broadens the absorption peaks. 

To calculate absorption-spectra we have to expand the system and add the ground 
state |0) to our basis. The Hamiltonian is thus extended to an 8 x 8 block matrix 

H'/ = I ? ^^ex 1 , (9) 






where we have chosen the zero exciton energy H^^ = and H^^^ is the 7x7 matrix 
defined in (|3]). 

Using the same block structure as in (jH]) the dipole matrix that governs the photon- 
exciton conversion is constructed as 

\ /^lox J 

where fi^^^ = fJ^i^^ and the zero to single excitation dipole block 

fe=l 

is defined using the seven dipole moments iik that give the excitation of BChlfc by the 
laser pulse. 

The absorption spectrum is calculated as 

/idM = // dte'-*tr(/i(t)MO)po)\ , (12) 

\^0 / rot 

where po = |0) (0| is the density operator for no excitation and we assume (5-shaped 
laser pulses. The different methods for performing the rotational average denoted by 
(■)rot are detailed below. The time evolution of the dipole operators in f|T2|) is given by 

In the FMO complex the seven BChl molecules have a fixed orientation with 
respect to each other. Within this configuration each BChl has a dipole vector d^ 
that characterises its interaction with light. We consider the dipole vectors of the 
BChls to be along the connection line of two nitrogen atoms Nb—^d [43j whose 
positions we obtain from the protein data bank [2^1 HI]. The norm of d^ accounts 
for the strength of the dipole and we assume d^ = d since all BChl molecules are 
identical. For the simulations we choose d = 1. Given a polarisation direction of the 
laser I = (sin /3 cos 7, sin /3 sin 7, cos where 7 G [0,27r) and /3 G [0,7r], the dipole 
moments in fill I) read 

fj,k = dk- 1 (13) 

Typical experiments are performed using spectroscopy of a sample solution and many 
FMO complexes are illuminated at the same time. Each of them is randomly orientated 
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with respect to the electric field in the laser pulse. If we assume a coherent dynamics 
expression f[T^ is readily evaluated in eigenbasis as 



(tr(/i(t)/x(0)po))rot = 




= E/^I.e-''^^*/'^ (15) 

where Ej is the jth eigenenergy and vl = {Ej\k) denotes the overlap of the kth state 
in site basis with the jth eigenstate. The integrationally averaged dipole moments in 
eigenbasis are 

/'27r P7T f \ 

'^^^ " ^ / / d/3 sin (3 i^dk- 1 vl) 
^ >^o Jo \k=l / 

= 3 XI ^i'^^^k ■ dm. (16) 

k,m=l 

Thus, we obtain the integrationally averaged squares of dipole moments in the eigenbasis 
/27rf2 = {0.25, 0.64, 0.39, 0.11, 0.62, 0.12, 0.20}, where the eigenenergies increase 
from left to right. For an appropriate d, these values are in rough agreement with 
the averaged dipole moments obtained from a fit to the experimental spectrum (i^yy — 
{49, 87, 73, 31, 82, 24, 36} debye^ [H]. 

However, we can also read the result as = | X^I m=i "^fc^m ^fc 'Udm- h, 
where the three polarisation vectors k are orthogonal to each other. This now implies 
that the rotational average can be performed using three orthogonal directions for I 
instead of sampling over the whole unit sphere. Note that this calculation is valid only 
if the interaction of system and bath does not interfere with the rotational average. 

In the following we discuss how the line-shapes of the peaks depend on the 
theoretical method and approximations commonly used for performing the time- 
dependent propagation in (|T2l) before we return to the different methods of performing 
the rotational average. 

3.1. Absorption spectra obtained by different propagation methods 

We calculate the absorption spectra of the FMO complex at the temperature T = 77 K. 
We propagate the exciton-bath system during 2048 fs with a numerical time step of 2 fs 
using the HEOM with reorganisation energy A = 35 cm^^ and a time scale for the bath 
correlations given by = 50 fs. Each HEOM propagation takes only 1 second on a 
GPU. 

In figure [T] (a), we compare the results obtained for the absorption spectrum 
calculated using the HEOM and the secular Born-Markov approximation (SBM) 
following Refs. The line shapes for each eigenenergy can be calculated in the 

eigenbasis Ij = e^'^* (tr [he^ (t) fJ-Ej (0) Po))^.^^- In the SBM simulation the complete 




Figure 1. (a) Absorption spectra /id(w) (fT^ for the FMO complex calculated at 
T — 77 K using the HEOM and the secular Born-Markov (SBM) approximation 
(without Lamb-shift) for 7"^ = 50 fs, A = 35 cm~^ and /ipiT ■ The vertical lines 
correspond to the eigenvalues of the exciton Hamiltonian ([3]). (b) Absorption lines of 
the second eigenstate calculated using the HEOM for different values of 7^^ in the 
spectral density 1^ with fixed A = 35 cm~^. 



spectrum is identical to the sum of the seven positive absorption peaks (not shown). Due 
to the rotating wave approximation the terms tr (/x^;^ (t) (0) Po) foi' 3 7^ ^ are zero. 
For the HEOM, these terms have non-vanishing contributions and the sum of the hne 
shapes do no longer yield the spectrum. Thus we have to calculate the line shape of 
the jth exciton peak as Ij = e^'^^ (ti (^^^i"" fJ'E, (t) ^Ej (0) Poj^ • Including these 
terms leads to more complicated line shapes with shoulders and negative contributions 
as shown in figure [1] (b) for the peak of the second eigenstate. We find that the more 
complicated line shapes arise from the non-secular approximations. Comparisons of 
different shapes of the spectrum obtained by different propagation methods and different 
spectral densities have been presented in j46[ Wl\ for the B850 ring of the light-harvesting 
complex LH2. 

We observe in figure [T] (a) that the energies for the HEOM peaks are shifted 
compared to the peaks obtained in the SBM approximation. Since we do not 
take into account the time-dependent Lamb shift in the SBM approximation (see 
eq. (8) in Ref. the corresponding SBM peaks are located at the exciton energies 
of Hamiltonian (|3]), marked by the vertical lines. Note that the eigenenergies of 
Hamiltonian ([3]) include the reorganisation energy A added to the "bare exciton- 
energies". Due to the dissipation of the reorganisation energy within the vibrational 
environment, taken into account by the HEOM [23], the peaks of the absorption 
spectrum shift to lower energies. For a fast phonon relaxation 7"^ — )■ (Markov limit) 
the reorganisation energy dissipates instantaneously and the peaks are shifted to lower 
energies by approximately A and then coincide with the bare exciton energies. Figure 
[T]^b) displays how the bath-correlation time affects the peak shapes and positions. For 
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Figure 2. (a) Absorption spectra /id(<^) of the FMO complex calculated using the 
HEOM with iVmax = 8 at T = 77 K, A = 35 cm'^ and = 166 fs. We compare 
different methods of performing the rotational average in the calculations. The dotted 
line results from using the fitted values /zfit, the dashed line takes the values j2. 
The solid line denotes a Monte Carlo (MC) average over 300 samples. The spectra 
are normalised with respect to the area under the graph, (b) Comparison of the 
absorption spectra for 7"^ = 50 fs for different values of static Gaussian disorder, which 
is characterised by its standard deviation (SD). In all cases the rotational average was 
performed with MC using 300 realisations. 



longer phonon relaxation times, the reorganisation energy dissipates slowly during the 
exciton dynamics and the peaks stay closer to the eigenenergies of ([3]). 

While the spectra shown above are calculated using fitted dipole moments from 
|14] . we will now analyse the changes due to rotational averaging on the absorption 
spectrum. In figure |2] (a) we show the results for the spectrum using the HEOM for 
A = 35 cm~^ and 7"^ = 166 fs. We compare different methods of rotational average 
discussed in the previous section, that is, a single propagation from integrationally 
averaged dipole moments using fi eq. ( IT6|) . a Monte Carlo (MC) set of 300 orientations 
and /ipiT- The results using three propagations with perpendicular laser directions 
xyz are not shown since the curve is almost identical to the MC spectrum for 300 
realisations. The MC method matches the experimental conditions closest, since the 
experiment shows the averaged spectrum of many randomly orientated complexes. To 
avoid propagating density matrices for many random orientations, the xyz approach 
can be chosen instead as both spectra show excellent agreement. Neither the integrated 
dipole moments, nor the fitted ones yield comparable results to the MC method and 
only in the Markovian limit it is possible to run the propagation with pre-averaged 
dipole moments. Here, we chose a longer bath correlation time 7"^ = 166 fs for the 
comparison of different rotational averages as the differences vanish for decreasing 7"^. 
In the following we always use several orientations for rotational averages and do not 
resort to pre-averaged dipole moments. 
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3.2. Broadening of the spectral lines 

If we compare the only rotationally averaged absorption spectra in figure [1] (a) with 
the experimentally measured ones [11] , we find that the rotationally averaged peaks are 
too narrow. Furthermore six peaks are well pronounced while in experiments several 
peaks overlap and the third exciton eigenstate is not individually resolved. This result 
is expected, since there are several additional effects besides rotational averaging, which 
broaden the peaks and bring the simulation results closer to the experimental data. If 
we increase either A (not shown) or the relaxation time 7"^, the peaks in the spectrum 
become broader, see figure [T](b). Both parameters change the spectral density and hence 
the coupling to the vibrational modes. The third mechanism is static disorder caused 
by the slowly fluctuating protein environment. We model the static disorder by adding 
a Gaussian distributed noise of a given standard deviation (SD) to each diagonal term 
in on]). The resulting spectra are shown in figure [2] (b) where we obtain broader peaks 
and show with increasing disorder the diminishing of the third peak in the spectrum. In 
addition we observe that the broadening changes the relative peak heights of the second 
and fourth peak yielding a better agreement with the experimental results. 

For the absorption spectrum changing the spectral density and inhomogeneous 
broadening yield similar results and it is not possible to distinguish which mechanism is 
relevant for the FMO complex. It is in the 2d echo-spectra where the difference becomes 
visible since the disorder results in a broadening along the diagonal frequencies only, 
while changing the exciton-phonon coupling broadens the peaks in all directions [18] . 



4. Two-dimensional echo-spectroscopy 

In addition to the information provided by absorption spectroscopy, two-dimensional 
echo-spectroscopy is used to study time-resolved processes such as energy transfer or 
vibrational decay as well as to measure intermolecular coupling strengths [IHl [501 EI] • 
Two-dimensional echo-spectroscopy consists of a sequence of three ultra-short pulses 
and a resulting signal pulse with fixed time-delay between the second and third pulse 
[52] [53] . The theory of 2d echo-spectroscopy is explained in detail in [53] [55] 156] . The 
2d echo-spectra are measured in two configurations which correspond to the rephasing 
(RP) and non-rephasing (NR) contributions of the third order response function 

S^'\t^, t2, tl) = e(t3)e(t2)e(ti) tr (/i (tg + t2 + tl) [/i (t2 + tl) , [/i (tl) , [/i(0), Po]]]) .(17) 

Using the impulsive approximation which assumes 5-functions for the temporal 
envelopes of the pulses, the two components of the spectrum corresponding to 
S^^Ht3,t2,ti) = 4?(t3,t2,ti) + 41(^3,^2,^1) read 

/OO POO 
dti / rft3e'"^*^-''^^*^4p(^3,t2,ti), (18) 
00 J —00 

/OO POO 
dti / rft3e'"^*^+'"^*^ 41(^3,^2,^1). (19) 
-00 J —00 
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Figure 3. The six Feynman diagrams included in the 2d echo-spectrum calculation. 
Vertical axis represent the time propagation and the action of the dipole operators are 
shown by diagonal arrows. For a complete theoretical description see [55] . 



Sorting this expression for the rephasing and the non-rephasing parts yields six different 
pathways, which represent the stimulated emission, the ground-state bleach and the 
excited state absorption [SU |55l EH |57], shown schematically in figure |3l By taking the 
time between second and third pulse as a control parameter as done in flTS]) and flTIJl) 
one can study the exciton dynamics which is probed by recording the emitted signal 
pulse. Due to the excited state absorption the dynamics is more involved compared to 
the one dimensional case. Including a double excited manifold, we construct an enlarged 
exciton Hamiltonian H^'^ which is block diagonal and of the form 

/ m' \ 



He 



2d 



(20) 




y i/|-y 

where H^^ = represents the ground state, H^'^^ is the single-exciton Hamiltonian with 
N^ites elements as defined in ([3]), and the energy levels of the two excitons are described 
by H^'^^. The two-exciton Hamiltonian H^'^ has [Ai"sites(-^sites ~ l)/2]^ entries which are 
constructed from the matrix elements of H^'^^ following [51} [H] . For the FMO complex 
H^'^ is a 29 X 29 matrix. In site basis, we denote the two exciton states by \ij), where 
< i < j < Ai'sites, and the diagonal entries of H^'^^ are given by 

{ij I I ^J) = {i I i/^ I i) + (j I if^ I j) . (21) 
The off-diagonal elements are given by 

I i/f ^ I kl) = (1 - {j I I /> + 6u (1 - 6,k) {j I i^e''" I k) 

+ (1 - Su) {i I Hl""^ I /) + 5ji (1 - 5,,) {i I Hl""^ I k) . (22) 
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Physically, this means that the \ij) and \kl) two-exciton states are coupled if they have 
one exciton in common while a second exciton is shifted from one site to another. The 
coupling constant is the same as for the transition of the corresponding exciton shift in 
the one exciton dynamics. Note also that the linear coupling to the independent baths, 
each of them associated to one BChl, leads to couplings between the one exciton state 
\i) and the two exciton states \ij) which do both influence the bath at site i and are in 
return influenced by the bath correlations. 

Following [31], we construct the dipole operator associated with the extended 
Hamiltonian (1201) 





















/^2ex 






(23) 



The zero to single excitation operators /i]';^ are defined as in ffTOj) and ffTTj) . The single to 
double excitation operators are defined according to (ij | /x^g^ | A;) = 6ik {j \ fJ.tex I O) + 

^jk {i I /^L I 0) and fi^^^ = {fitexV . 

In order to calculate 2d spectra, we have to evolve the 29 x 29 density matrix to 
obtain the third order response function S^^\ The hierarchical method incorporates the 
effect of the baths through the auxiliary system of the a-matrices. When calculating 
the 2d spectra, we propagate according to the Feynman diagrams shown in figure [3] and 
act on the density matrix and the auxiliary matrices with dipole operator to simulate 
the interaction with the laser field. 

We observe in (|T7|) that the 2d spectra depend on the fourth order dipole moments. 
To take into account the random orientation of the samples we need to compute 
rotational averages {lAl^'f)^ot ^^^^ performed by sampling 10 laser polarisation- 
vectors aligned to the vertices of a dodecahedron in one half of the coordinate space. 
We take the vertices (±1,±1,1), ^O,±i,0j, ^±^,0,0^ and ^±0,0,^^, where = 

I (l + y/h) denotes the golden ratio. We have compared the dodecahedral summation 
with the full spherical averaging using a 10000-shoot MC and obtain differences of less 
than 1 percent of the MC results. 



Numerical results for the 2d echo-spectrum of the FMO complex 

For the simulation of 2d spectra, we calculate the third order response function 
S^^\tiit2-,t'i) in (fT7|) . We obtain 128^ data points with a spacing of 16 fs in ti and 
direction using a numerical time step of 2 fs. This means that we have to perform 
over 131 000 propagation steps for each of the six Feynman diagrams. Depending on 
the delay time T^eiay = ^2 our simulations take 48 to 53 minutes on a single NVIDIA 
Fermi C2050 GPU which provides a 38 fold speed up compared to the CPU algorithm 
with filtering methods that reduce the hierarchy ^17]. We always show the sum of the 
real parts of the rephasing and the non-rephasing spectrum. 

Figure H] (a) shows a simulated 2d spectrum for a rotationally averaged sample 
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Figure 4. (a) Rotationally averaged 2d spectrum i?e[/NR + /Rp] of the FMO complex 
calculated using the HEOM with A = 35 cm"\ T = 77 K and = 50 fs and a 
delay time Tdeiay — fs. (b) Oscillations as a function of the delay time for the (1, 1) 
and (5, 5) diagonal peak and the (1, 5) peak below the diagonal. The points have been 
connected to guide the eye. 

with no disorder at T = 77 K. The figure is encoded using a hnear colour scale. We 
divided the raw data by a constant to obtain an amplitude on a —10-40 range. The 
thin white lines have been added to indicate the expected peak positions. The peaks 
are denoted by their excitonic energy-level indices {i,j), with the first number referring 
to the cj3-axis and the second one to the wi-axis. The spectrum shows similarities to 
the experimentally measured one in [H] insofar as it reproduces the large peak on the 
diagonal at the second eigenstate and a second near the fourth and fifth. Furthermore 
we find negative amplitudes above the diagonal and positive peaks below. As expected 
from the discussion of the absorption spectrum in section 13. 2[ the peak shapes differ 
from the experimental ones due to the lack of static disorder. 

Figure H] (b) shows the peak oscillations for the diagonal peaks (1,1) and (5,5) 
and the off-diagonal peak (1,5) below the diagonal. We renormalised the peak heights 
at Tdeiay = fs to Unity such that the time resolved amplitude of the three peaks are 
conveniently put into the same panel. We choose these peaks because figure H] shows 
peak (1, 5) as the highest off-diagonal peak in the line where U3 equals the lowest lying 
eigenstate. We see that diagonal and off-diagonal peaks oscillate in time. This indicates 
that our model imitates the coherent energy transport within the FMO complex that 
has been observed experimentally O [6]. Theoretically, an off-diagonal peak {i,j) is 
expected to undergo oscillations with a frequency corresponding to the difference of its 
excitonic eigenenergies {Ej — Ei). For the (1,5) peak we expect a period of 99 fs which 
is in good agreement to our numerical result yielding a best-fit period of 97 fs. When 
we calculate the peak heights, we integrate over squares of width Aw = 16 cm~^ centred 
around the white crosses. We follow the changes of the peak heights with increasing 
delay time by plotting the absolute value of the spectrum normalised with respect to 



Modelling of Oscillations in 2d Echo-Spectra of the Fenna- Matthews- Olson Complex 15 



(a) (b) 




Figure 5. (a) Amplitude of tlie cross peak (1, 5) of tlie rotationally averaged spectra of 
tlie FMO coniplex as a function of delay time for different temperatures. Temperature 
dependence of the peak oscillations for the (1,5) cross peak, (b) Oscillations in the 
peak width of the (1,1) diagonal peak measured in diagonal and anti-diagonal direction 
compared to the peak height at T = 77 K. The points have been connected to guide 
the eye. 

the value at Tdeiay = fs. 

In order to analyse the robustness of the peak oscillations when the thermal 
fluctuations increase we plot in flgure |5] (a) the oscillations of the (1,5) cross peak 
for several temperatures. We compare with the experimental data in [6] and observe a 
qualitative similar behaviour with weaker oscillations and faster decay when temperature 
is increased. We do not obtain the picoseconds lasting beatings that have been 
observed experimentally. One possible factor is the shape of the spectral density that 
is used for the HEOM. We flt the oscillations in flgure [5] (a) by an exponentially 
decaying sine. The flt is done separately for the stimulated emission and excited 
state absorption pathway, which couple to one or two baths respectively. Based on 
the damping of the oscillations on peak (1,5) the temperature dependence of the 
dephasing rate of our simulation is given by 7feph iT)/T = (0.66 ± 0.05) cm~^/K and 
7d^h (^)/^ = (0.50 ± 0.03) cm~^/K. This is in good agreement to what we expect for 
the chosen spectral density (2/csA/(7T) = 0.46 cm~^/K) and the value 0.52 cm~^/K 
determined from the experimentally measured data [6]. 

In the experiments an ant i- correlation in peak shapes between the diagonal peak 
width and the height |5J, or alternatively the ratio of the diagonal peak- width over the 
anti-diagonal height [4J has been found. The width was measured by fltting a Gaussian 
to a diagonal cut through the spectrum and taking its standard deviation. Applying the 
same procedure to the simulation data yields also anti-correlations shown in flgure |5](b). 

Next, we analyse the role of static disorder in the results for the 2d echo-spectra. 
When including static disorder, the peaks get elongated along the diagonal as observed 
when comparing the spectra at zero delay time shown in flgures S] (a) and [6] (a) and (b) 
for increasing values of disorder. The plots with disorder show better agreement with 
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Figure 6. Two-dimensional echo-spectra i?e[/NR, + -^rp] of the FMO complex with 
rotational and disorder averages at temperature T = 77 K. The averages are performed 
over 400 realisations. The spectral-density parameters are A = 35 cm~^ and = 
50 fs. (a) Gaussian disorder with standard deviation (SD) 25 cm~^ at delay time 

Tdelay = fs. (b) SD= 50 Cm^^, Tdelay = fs. (c) SD= 25 Cm"!, Tdclay = 4000 fs. 

the experiment [H]. Figure [6] (c) shows the 2d spectra for SD= 25 cm~^ for a long delay 
time Tdelay = 4000 fs for which the system approaches its thermal stationary state. We 
observe that the excitation moves from the diagonal peaks into the regions with lower 
frequency u^. This occurs as the system relaxes to lower lying eigenstates during the 
delay time. In experiments by Brixner et al [9] the decay of the population into the 
ground state of the system is also visible. 

The absence of the experimentally seen oscillation in the HEOM calculation 
in Ref. [17] has led Chen et al to conclude that the importance of static disorder 
requires careful reconsideration. First-principles simulations of static disorder in the 
FMO complex are computationally extremely time-consuming, since several hundred 
of realisations are required to reach converged results. In order to study the role of 
static disorder in the oscillations of the peaks of the spectrum and to ensure numerical 
convergence, we construct a dimer to model the (1, 5) peak of the FMO spectrum. For 
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Figure 7. Peak oscillations in a model dimer for T = 77 K, A = 35 cm^^ and 
= 50 fs. To concentrate on the disorder effects, no rotational average is performed 
and the calculation is done using = —5.4 and ^2 — 12.2 (a) no static disorder, (b) 
Gaussian disorder with standard deviation (SD) 25 cm~^. The average was done over 
1200 realisations at each delay time. The points have been connected to guide the eye. 



the model dimer we chose the Hamiltonian such that the first and fifth FMO complex 
eigenvalues are reproduced. The couphng between both sites is set to 50 cm^^. The 
simulations with and without disorder show oscillations on the diagonal and off-diagonal 
peaks. The most striking effect of disorder is a decreased oscillatory amplitude. While 
the cross peak in the simulation without disorder has lost 45 percent of its amplitude in 
the first minimum, in the disorder case the loss is only 35 percent. From a fit, we observe 
that the periods of the disorder curves are slightly shorter T = 92 fs and T = 98 fs for 
the (1, 1) and (2, 1) curves while in the no-disorder plot T = 95 fs and T = 98 fs has been 
obtained, respectively. For the FMO complex we expect that oscillations are reduced in 
amplitude by a similar magnitude and therefore still persist within the model. 

5. Conclusions 

We have presented calculations of the absorption spectra and the 2d echo-spectra for the 
FMO complex at different temperatures using the HEOM for solving the exciton and 
the vibrational bath dynamics and adding static disorder to account for the fiuctuating 
environment. Important ingredients necessary to find good agreement with experiments 
for the absorption spectra using the HEOM propagation method is to perform the 
rotational average and to include static disorder. We have also discussed how other 
simplified propagation methods affect the line shapes and shift the peaks in the spectra. 

Our rotationally averaged 2d echo-spectra reflect the main features of the 
experimentally measured ones. At delay time T^eiay = we flnd negative contributions 
above the diagonal and elongated peaks in the diagonal which are transferred to the 
lowest exciton states after long times on the picosecond scale. The HEOM method 
yields an oscillatory behaviour of peak amplitudes over a wide temperature range. At 
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T = 77 K we obtain oscillations of cross peak amplitudes up to 300 fs with the expected 
frequency that corresponds to the difference in excitonic eigenenergies. In addition 
oscillations on the diagonal peaks are present. The long duration of oscillations in 
the experimental data (up to 1.8 ps in Ref. [B]) suggests that further refinements of 
the theoretical models are necessary, for example by adjusting the form of the spectral 
density which affects the population dynamics |38] . 

Remarkably, we find excellent agreement with the measured damping rate of the 
oscillations due to thermal fluctuations. This implies a good choice of the parameters 
of the model and might explain why we obtain different results compared to previous 
calculations, which were done for a higher value of = 100 fs [T7], leading to an 
increased dephasing by a factor of two. At ambient temperatures the HEOM have the 
advantage to require less terms to converge while approximate methods such as the 
secular and full Redfield approximations do not yield the correct time-evolution of the 
density matrix and overestimate the thermalisation rate. 

In agreement with the experimental measurements we observe the anti-correlated 
oscillations of the peak heights and peak shapes for the diagonal peaks. The analysis 
of a model dimer shows that the cross peak oscillations persist even in the presence of 
static disorder which mainly results in a reduction of oscillation amplitude by about 30 
percent. 
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